R Notebook

Exploring dual-cue approach

library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(ggplot2)
library(here)

here() starts at /Users/avrilwang/Desktop/Project-Plasmodium

library(deSolve)
library(crone)
library(optimParallel)

Loading required package: parallel

library(doParallel)

Loading required package: foreach
Loading required package: iterators

library(doRNG)

Loading required package: rngtools

library(arrow)

Attaching package: ‘arrow’

The following object is masked from ‘package:utils’:

    timestamp

library(stringr)
library(parallel)
library(ggpubr)
library(microbenchmark)
library(profvis)

Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

library(RANN)

parameters_tsukushi <- c(R1 = 8.89*(10^6), # slightly higher
                lambda = 3.7*(10^5),
                mu = 0.025, 
                p = 8*(10^-6), # doubled form original
                alpha = 1, 
                alphag = 2, 
                beta = 5.721, 
                mum = 48, 
                mug = 4, 
                I0 = 43.85965, 
                Ig0 = 0, 
                a = 150, 
                b = 100, 
                sp = 1,
                psin = 16.69234,
                psiw = 0.8431785,
                phin = 0.03520591, 
                phiw = 550.842,
                iota = 2.18*(10^6),
                rho = 0.2627156)

I_range_none <- seq(0, 6*(10^6), by = (6*(10^6))/500)
I_range_log <- seq(0, log10(6*(10^6)), by = (log10(6*(10^6)))/500)
G_range_none <- seq(0, 6*(10^4), by = (6*(10^4))/500)
G_range_log <- seq(0, log10(6*(10^4)), by = (log10(6*(10^4)))/500)

testing out mgcv dual cue spline surface. looks good.

note that if we go with ti, than we have only interactions between the 2 cues whereas te includes both main effect and interaction. ti gives us 5 parameters to optimze whereas te gives us 9. We should compare the two

source(here("functions/chabaudi_si_clean.R"))

# smaller division
cue_range <- seq(0, 6*(10^6), by = (6*(10^6))/500)
cue_range_b <- seq(0, 6*(10^4), by = (6*(10^4))/500)
parameters_cr <- rep(-0.5,9)

cr_grid <- expand.grid(cue_range, cue_range_b)
    ## rename
    names(cr_grid) <- c("cue_range", "cue_range_b")
    ## create dummy y
    dummy_y <- runif(length(cue_range_b), 0, 1)
    ## put together df
    dummy_df <- data.frame(cue_range, cue_range_b, dummy_y)
    ## gam model
    dummy_cr.mod <- mgcv::gam(dummy_y ~ te(cue_range, cue_range_b, 
                                           bs = c("tp", "tp")
                                           k = c(3,3)), 
                              data = dummy_df, 
                              method = "REML")
    ## assign parameters
    dummy_cr.mod$coefficients <- parameters_cr
    
    # exponential transformation to limit conversion rate to between 0 and 1
    cr_fun <- function(cue_1, cue_2){
      res <- exp(-exp(mgcv::predict.gam(dummy_cr.mod, 
                                        newdata = data.frame("cue_range" = cue_1,
                                                             "cue_range_b" = cue_2))))
      return(res)
    }
    

cutting down computational time

source(here("functions/chabaudi_si_clean.R"))
source(here("functions/test.R"))

# each run takes around 1 minute (cr_fun) (100 divisions)
# switching to nn2 takes 145 seconds
# increasing divisons to 500 using predict.gam takes still around 57 seconds. So increasing to 500 doesn't seem like a bad idea
system.time(test(parameters_cr = rep(0.5, 5),
       immunity = "tsukushi",
       parameters = parameters_tsukushi,
       time_range = seq(0, 20, 0.001),
       cue = "I",
       cue_b = "G",
       cue_range = I_range_log,
       cue_range_b = G_range_log,
       solver = "vode",
       log_cue = "log10",
       log_cue_b = "log10"
       ))

making shit go fast! Need to cut down on time to get cr from function

cue_range <- I_range_log
cue_range_b <- G_range_log
cr_grid <- expand.grid(cue_range, cue_range_b)

## rename
names(cr_grid) <- c("cue_range", "cue_range_b")
## create dummy y
dummy_y <- runif(length(cue_range_b), 0, 1)
## put together df
dummy_df <- data.frame(cue_range, cue_range_b, dummy_y)
## gam model
dummy_cr.mod <- mgcv::bam(dummy_y ~ ti(cue_range, cue_range_b, 
                                           k = c(3,3)), 
                              data = dummy_df)
## assign parameters
dummy_cr.mod$coefficients <- rep(0.5,5)
    
# original function of getting cr
cr_fun <- function(cue_1, cue_2){
      exp(-exp(mgcv::predict.bam(dummy_cr.mod, 
                                 newdata = data.frame("cue_range" = cue_1,
                                                             "cue_range_b" = cue_2),
                                 n.threads = 3)))
  }
    

cr_res <- exp(-exp(mgcv::predict.gam(dummy_cr.mod, 
                                        newdata = cr_grid)))
cr_df <- cbind(cr_grid, cr_res)

cue_1 <- 1
cue_2 <- 0.5

# another function for nearnest neighoubrs
nn_fun <- function(cue_1, cue_2){
  cr_df[nn2(data=cr_grid, query = cbind(cue_1, cue_2), k=1, treetype= "kd",
            searchtype = "priority")[[1]],3]
}

nn_fun2 <- function(cue_1, cue_2){
  cr_df[FNN::knnx.index(data=cr_grid, query = cbind(cue_1, cue_2), k = 1, algorithm = "kd_tree"),3]
}


# very similar
cr_fun(1,1)
nn_fun(1,1)
nn_fun2(1,1)

# gam.predict is faster
microbenchmark(cr_fun(1,1))
microbenchmark(nn_fun(1,1))
microbenchmark(nn_fun2(1,1))

sanity checkz

source(here("functions/test.R"))
# switching between cue and cue_b. should not make a difference
# 4.469067
test(parameters_cr = rep(0.5, 5),
       immunity = "tsukushi",
       parameters = parameters_tsukushi,
       time_range = seq(0, 20, 0.01),
       cue = "I",
       cue_b = "G",
       cue_range = I_range_log,
       cue_range_b = G_range_log,
       solver = "vode",
       log_cue = "log10",
       log_cue_b = "log10"
       )

# 4.469066. good
test(parameters_cr = rep(0.5, 5),
       immunity = "tsukushi",
       parameters = parameters_tsukushi,
       time_range = seq(0, 20, 0.01),
       cue = "G",
       cue_b = "I",
       cue_range = G_range_log,
       cue_range_b = I_range_log,
       solver = "vode",
       log_cue = "log10",
       log_cue_b = "log10"
       )

# short time steps vs longer time steps. barely made a difference
# 4.46904
test(parameters_cr = rep(0.5, 5),
       immunity = "tsukushi",
       parameters = parameters_tsukushi,
       time_range = seq(0, 20, 0.001),
       cue = "I",
       cue_b = "G",
       cue_range = I_range_log,
       cue_range_b = G_range_log,
       solver = "vode",
       log_cue = "log10",
       log_cue_b = "log10"
       )

# changing logs on cue_b
# 6.20639
test(parameters_cr = rep(0.5, 5),
       immunity = "tsukushi",
       parameters = parameters_tsukushi,
       time_range = seq(0, 20, 0.01),
       cue = "I",
       cue_b = "G",
       cue_range = I_range_log,
       cue_range_b = G_range_none,
       solver = "vode",
       log_cue = "log10",
       log_cue_b = "none"
       )

# changing logs on cue
# 5.110873
test(parameters_cr = rep(0.5, 5),
       immunity = "tsukushi",
       parameters = parameters_tsukushi,
       time_range = seq(0, 20, 0.01),
       cue = "I",
       cue_b = "G",
       cue_range = I_range_none,
       cue_range_b = G_range_log,
       solver = "vode",
       log_cue = "none",
       log_cue_b = "log10"
       )

# checing single cue infection. same 
## 4.724804
test(parameters_cr = rep(0.5, 4),
       immunity = "tsukushi",
       parameters = parameters_tsukushi,
       time_range = seq(0, 20, 0.01),
       cue = "I",
       cue_range = I_range_none,
       solver = "vode",
       log_cue = "none",
       )
# 4.724804
chabaudi_si_clean(parameters_cr = rep(0.5, 4),
       immunity = "tsukushi",
       parameters = parameters_tsukushi,
       time_range = seq(0, 20, 0.01),
       cue = "I",
       cue_range = I_range_none,
       solver = "vode",
       log_cue = "none",
       )

#————————————————–# # remaking of co-infection reaction norm #————————————————–#

source(here("functions/chabaudi_ci_clean.R"))
source(here("functions/par_to_df.R"))

function for getting conversion rate reaction norm and dynamics

get_ci_dyn <- function(df){
  
  # get heavisde transformation
  heaviside_trans <- function(cue_range, max){
    crone::heaviside(cue_range)*(cue_range)+(crone::heaviside(cue_range-max)*(max-cue_range))
}
  
  # read in df info
  cue_1 <- ifelse(stringr::str_detect(df$cue, "-i"), 
                  gsub("*-i", "1", df$cue),
                  df$cue)
  
  cue_2 <- ifelse(stringr::str_detect(df$cue, "-i"), 
                  gsub("*-i", "2", df$cue),
                  df$cue)
  
  log <- ifelse(stringr::str_detect(df$log, "log"),
                "log10", "none")
  
  # get parameter
  par <- c(df$var1, df$var2, df$var3, df$var4)
  cue_range <- seq(df$low, df$high, by = df$by)
  label <- df$label
  
  # get dynamics
  dyn <- chabaudi_ci_clean(parameters_cr_1 = par,
                  parameters_cr_2 = par,
                  immunity = "tsukushi",
                  parameters = parameters_tsukushi,
                  cue_1 = cue_1,
                  cue_2 = cue_2,
                  cue_range_1 = cue_range,
                  cue_range_2 = cue_range,
                  log_cue_1 = log,
                  log_cue_2 = log,
                  solver = "vode",
                  time_range = seq(0, 30, 0.001),
                  dyn = T)
  
    dyn2 <- data.frame(dyn, label = rep(label, nrow(dyn)))

  # get reaction norm
  rn <- par_to_df(par = par, cue_range = cue_range)

  # get rug
  ## if + in cue, need to filter out those variables and add together
  if(stringr::str_detect(cue_1, "\\+")){
    cue_split <- stringr::str_split(string = cue_1, pattern = "\\+", simplify = T)
    ## get the two cues 
    cue_temp_1 <- cue_split[[1]]
    cue_temp_2 <- cue_split[[2]]
    
    ## filter dyn
    rug <- dyn2 %>% 
      filter(variable == cue_temp_1 | variable == cue_temp_2) %>% 
      dplyr::group_by(time) %>% 
      dplyr::mutate(sum = sum(value)) %>% 
      select(time, value = sum, label)
  }
  
  ## if sum
  if(cue_1 == "sum"){
    ## filter dyn
    rug <- dyn2 %>% 
      filter(variable == "I1" | variable == "I2" | variable == "Ig1" | variable == "Ig2") %>% 
      dplyr::group_by(time) %>% 
      dplyr::mutate(sum = sum(value)) %>% 
      select(time, value = sum, label)
  }
  
  # for simplified cue
  if(stringr::str_detect(cue_1, "\\+", negate = T) && cue_1 != "sum"){
    rug <- dyn2 %>% 
      dplyr::filter(variable == cue_1) %>% 
      dplyr::select(time, value, label)}
  
  # for logged, backtransform cue range of reaction norm 
  rn2 <- rn
  
  if(stringr::str_detect(log, "log")){rn2$cue_range <- 10^(rn2$cue_range)}

  # append label to reaction norm, dyn, and rug
  rn2 <- data.frame(rn2, label = rep(label, nrow(rn)))

  return(list(dyn2, rn2, rug))
}

read in co_infection opt

ci_opt.df <- read.csv(here("data/ci_opt.csv"))
cue_range_ci <- read.csv(here("data/cue_range_ci.csv"))

# join with cue range
ci_opt.df <- ci_opt.df %>% left_join(select(cue_range_ci, id, low, high, by), by = c("id" = "id"))

# split
ci_opt.ls <- split(ci_opt.df, seq(nrow(ci_opt.df)))

# loop
ci_opt.dyn <- lapply(ci_opt.ls, get_ci_dyn)

# bind together dynamics
dyn.ls <- lapply(ci_opt.dyn, function(x)x[[1]])
ci_dyn.df <- do.call(rbind, dyn.ls)
write_parquet(ci_dyn.df, here("data/ci_dyn/ci_dyn.parquet"))

rn.ls <- lapply(ci_opt.dyn, function(x)x[[2]])
ci_rn.df <- do.call(rbind, rn.ls)
write_parquet(ci_rn.df, here("data/ci_dyn/ci_rn.parquet"))

rug.ls <- lapply(ci_opt.dyn, function(x)x[[3]])
ci_rug.df <- do.call(rbind, rug.ls)
write_parquet(ci_rug.df , here("data/ci_dyn/ci_rug.parquet"))

get in single infection reaction norm to match

si_opt.df <- read.csv(here("data/si_opt.csv"))
cue_range_si <- read.csv(here("data/cue_range_si.csv"))

si_opt.df <- si_opt.df %>% left_join(select(cue_range_si, id, low, high, by, label), by = "id")

# function to get single infection reaction norm
get_si_rn <- function(df){
  par <- c(df$var1, df$var2, df$var3, df$var4)
  cue_range <- seq(df$low, df$high, by = df$by)
  
  rn <- par_to_df(par = par, cue_range = cue_range)
  
  # for logged, backtransform cue range of reaction norm 
  rn2 <- rn
  if(stringr::str_detect(df$log, "log")){rn2$cue_range <- 10^(rn2$cue_range)}

  # append label to reaction norm, dyn, and rug
  rn2 <- data.frame(rn2, id = df$id, label = rep(df$label, nrow(rn)))
  
  return(rn2)
}
# rn function and bind together
si_opt.ls <- split(si_opt.df, seq(nrow(si_opt.df)))
si_opt.rn <- lapply(si_opt.ls, get_si_rn)
si_opt.rn2 <- do.call(rbind, si_opt.rn)

get single infeciton rug

si_dyn <- read_parquet(here("data/si_dyn/si_dyn_30.parquet"))

# join to get label
si_dyn.df <- left_join(si_dyn, si_opt.df, by = "id")

# process get rug
get_rug <- function(df){
  cue <- unique(df$cue)
  
   if(stringr::str_detect(cue, "\\+")){
    cue_split <- stringr::str_split(string = cue, pattern = "\\+", simplify = T)
    ## get the two cues 
    cue_temp_1 <- cue_split[[1]]
    cue_temp_2 <- cue_split[[2]]
    ## filter dyn
    rug <- df %>% 
      filter(variable == cue_temp_1 | variable == cue_temp_2) %>% 
      dplyr::group_by(time) %>% 
      dplyr::mutate(sum = sum(value)) %>% 
      select(time, value = sum, label)
  }
  
  # for simplified cue
  if(stringr::str_detect(cue, "\\+", negate = T)){
    rug <- df %>% 
      dplyr::filter(variable == cue) %>% 
      dplyr::select(time, value, label)}
  
  return(rug)
}

# split based on individual labels
si_dyn.ls <- si_dyn.df %>% group_split(id)

# run function to get rug
si_rug <- mclapply(si_dyn.ls, get_rug)

si_rug.df <- do.call(rbind, si_rug)

process reaction norm data to limit reaction norm to around the infection space

match single infection rn with coinfection

plot reaction norm

c

function (...)  .Primitive("c")

#—————————————–# # plot maximum fitness accumulation for both si and ci #—————————————–#

ez_label <- read.csv(here("data/ez_label.csv"))

process single infection data

# import in si maximum transmission potential
si_dyn.df <- read_parquet(here("data/si_dyn/si_dyn_30.parquet"))

# get maximum tau_cum for 30 days
si_max_tau <- si_dyn.df %>% 
  filter(variable == "tau_cum") %>% 
  group_by(id) %>% 
  summarise(fitness_si = max(value))

process co-infection data

# import in co-infeciton data, 30 days simulation
ci_dyn.df <- read_parquet(here("data/ci_dyn/ci_dyn.parquet"))

# get maximum tau_cum. Note multiplied by 2 to get total fitness for both strains, given that both are sadiptng the same strategy, we consider them the same genotye
ci_max_tau <- ci_dyn.df %>% 
  filter(variable == "tau_cum1") %>% 
  group_by(label) %>% 
  summarise(fitness_ci = max(value, na.rm = T))

get label and join together 2 dataframes

si_ci_max_tau <- ci_max_tau %>% 
  left_join(ez_label, by = c("label" = "label_ci")) %>% 
  left_join(si_max_tau, by = c("id_si" = "id"))

plot

ggsave(here("figures/report21/si_ci_fitness.png"), height = 5)

Saving 7 x 5 in image

#—————————————# # time series conversion rate (line form) #—————————————#

single infection

cr_si <- si_dyn.df %>% 
  filter(variable == "cr")

# good cue bad cue
si_cue_dv <- si_max_tau %>% 
  mutate(classification = case_when(
    fitness_si > 10.5 ~ "High-performing",
    fitness_si <= 10.5 ~ "Poor-performing"
  ))

# join with classificaiton
cr_si.df <- cr_si %>% left_join(si_cue_dv) %>% left_join(ez_label, by = c("id" = "id_si"))

Joining, by = "id"

# split into top erforming and poor-performing cues
cr_si.high <- cr_si.df %>% filter(classification == "High-performing")
cr_si.poor <- cr_si.df %>% filter(classification == "Poor-performing")

# plot
cr_si_pl.pre <- ggplot() +
  geom_line(data = cr_si.poor, aes(x= time, y = value, group = id), color = "dark grey") +
  labs(color = "Single infection\nhigh-performing cues", x = "Time (days)", y = "Conversion rate") +
  theme_bw()+ scale_y_continuous(labels = function(x) format(x, scientific = TRUE, accuracy = 0.1)) +
  guides(shape = FALSE)

Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.

cr_si_pl <- cr_si_pl.pre +
  geom_line(data = cr_si.high, aes(x= time, y = value, group = ez_label_si, color = ez_label_si), size = 1) +
  scale_color_manual(values=c( "#4575b4", "#fc8d59", "#fdcb44", "#91bfdb"))

co-infection

# get relevent variables
cr_ci <- ci_dyn.df %>% 
  filter(variable == "cr_1")

# good cue bad cue
ci_cue_dv <- ci_max_tau %>% 
  mutate(classification = case_when(
    fitness_ci > 2.25 ~ "High-performing",
    fitness_ci <= 2.25 ~ "Poor-performing"
  ))

si_cue_dv

# join with classificaiton
cr_ci.df <- cr_ci %>% left_join(ci_cue_dv )%>% left_join(ez_label, by = c("label" = "label_ci"))

# split into top erforming and poor-performing cues
cr_ci.high <- cr_ci.df %>% filter(classification == "High-performing")
cr_ci.poor <- cr_ci.df %>% filter(classification == "Poor-performing")

# plot
cr_ci_pl.pre <- ggplot() +
  geom_line(data = cr_ci.poor, aes(x= time, y = value, group = label), color = "dark grey", alpha = 0.5) +
  labs(color = "Co-infection\nhigh-performing cues", x = "Time (day)", y = "Conversion rate") +
  theme_bw()+ scale_y_continuous(labels = function(x) format(x, scientific = TRUE, accuracy = 0.1)) +
  guides(shape = FALSE)

cr_ci_pl <- cr_ci_pl.pre +
  geom_line(data = cr_ci.high, aes(x= time, y = value, group = ez_label, color = ez_label), size = 1) +
  scale_color_manual(values=c( "#4575b4", "#fc8d59", "#fdcb44", "#91bfdb"))

arrange together

cr_pl <- ggarrange(cr_si_pl, cr_ci_pl, align = "hv", ncol = 1)

#—————————————# # Time series conversion rate (heatmap) #—————————————# note that the order of cues is sorted by descending 30-days fitness. # single infection

cr_hm_si.pl1 <- ggplot() +
  geom_tile(data = cr_si.poor, aes(x = time, y = forcats::fct_reorder(ez_label_si, fitness_si), fill = value)) +
  labs(x = "Time (days)", y = "Low performing\nsingle infection cues", fill = "Conversion rate") +
  scale_fill_viridis_c(limits = c(0, 1)) +
  xlim(1, 30) +
  theme_bw()

cr_hm_si.pl2 <- ggplot() +
  geom_tile(data = cr_si.high, aes(x = time, y = forcats::fct_reorder(ez_label_si, fitness_si), fill = value)) +
  labs(x = "", y = "High performing\nsingle infection cues", fill = "Conversion rate") +
  scale_fill_viridis_c(limits = c(0, 1)) +
  xlim(1, 30) +
  theme_bw()

co-infection

cr_hm_ci.pl1 <- ggplot() +
  geom_tile(data = cr_ci.poor, aes(x = time, y = forcats::fct_reorder(ez_label, fitness_ci), fill = value)) +
  labs(x = "Time (days)", y = "Low performing\nco-infection infection cues", fill = "Conversion rate") +
  scale_fill_viridis_c(limits = c(0, 1)) +
  xlim(1, 30) +
  theme_bw() +
  theme(legend.position = "none")

cr_hm_ci.pl2 <- ggplot() +
  geom_tile(data = cr_ci.high, aes(x = time, forcats::fct_reorder(ez_label, fitness_ci), fill = value)) +
  labs(x = "", y = "High performing\nco-infection cues", fill = "Conversion rate") +
  scale_fill_viridis_c(limits = c(0, 1)) +
  xlim(1, 30) +
  theme_bw() +
  theme(legend.position = "none")


cr_hm.pl <- ggarrange(cr_hm_si.pl2, cr_hm_ci.pl2, cr_hm_si.pl1, cr_hm_ci.pl1, align = "hv", ncol = 2, nrow = 2, common.legend = T)

Warning: Removed 6000 rows containing missing values (geom_tile).
Warning: Removed 6000 rows containing missing values (geom_tile).
Warning: Removed 4000 rows containing missing values (geom_tile).
Warning: Removed 10000 rows containing missing values (geom_tile).
Warning: Removed 12000 rows containing missing values (geom_tile).

ggsave(here("figures/report21/conversion_rate_heatmap.png"), height = 5.5, width = 12)

#—————————————# # disase curve #—————————————# note we are using the same 30 days infection dynamics from before

for single infection

note that we are dividing the cues into good and bad based on arbitrary cut off at 10.5

co-infection

# plot both single and co-infeciton disease curve together

disease_curve_pl <- ggarrange(disease_curve_pl.1, disease_curve_pl.2, align = "hv", heights = c(0.75, 1))

Error in ggarrange(disease_curve_pl.1, disease_curve_pl.2, align = "hv",  : 
  object 'disease_curve_pl.1' not found

#———————–# Arranging opt_cue figure #———————–#

ggarrange(opt_cue_pl.A, opt_cue_pl.BC, labels = c("A", ""), ncol = 2, widths = c(0.75, 1))

Warning: ggrepel: 7 unlabeled data points (too many overlaps). Consider increasing max.overlaps

ggarrange(opt_cue_pl.A, opt_cue_pl.BC, labels = c("A", ""), ncol = 2, widths = c(0.75, 1))

Warning: ggrepel: 7 unlabeled data points (too many overlaps). Consider increasing max.overlaps

ggsave(here("figures/report21/opt_cue_final.png"), width = 14, height = 9)

#———————————# # remaking heatmaps #———————————# #———————————-# # static heatmap #———————————-# # processing data

max(static.df4$fitness_diff)

[1] 0.9096023

get heatmap

ggsave(here("figures/report21/static_heatmap.png"), width = 10)

Saving 10 x 7 in image

#————————-# # ci invasion matrix #————————-#

# import and process ci-invasion data
invade.ls <- list.files(here("data/ci_invasion/"), pattern = "*.csv", full.names = T)

invade.ls2 <- lapply(invade.ls, read.csv)
invade.df <- do.call(rbind, invade.ls2)
#write.csv(invade.df, row.names = F, here("data/ci_invasion.csv"))
invade.df <- read.csv(here("data/ci_invasion.csv"))

sanity check. all invasion combination hsould have a fitness higher than static competition

sc_static.df <- static.df %>% 
  mutate(fitness_static = fitness_1 - fitness_2) %>% 
  select(id_1, id_2, fitness_static)

# get reverse
sc_static.df2 <- sc_static.df
sc_static.df2$fitness_static <- sc_static.df$fitness_static*-1
names(sc_static.df2) <- c("id_2", "id_1", "fitness_static")

# combine
sc_static.df3 <- rbind(sc_static.df, sc_static.df2)

# left join. none of these are huge and invasion runs for 20 days whereas static competition runs for 30 days. looks good
invade.df %>% 
  left_join(sc_static.df3, by = c("mut_id"= "id_1", "res_id" = "id_2")) %>% 
  mutate(invade_static = fitness - fitness_static) %>% 
  filter(invade_static < 0)

process data for invasion matrix

create summary bar chart

plot togehter

invasion.pl <- ggarrange(invasion.pl1, invasion.pl2, align = "h", common.legend = T, widths = c(2, 1))
ggsave(here("figures/report21/invasion_heatmap.png"))

Saving 7.29 x 4.51 in image

#————————————-# # Effect of various cue perception #————————————-#

#——————————-# # static #——————————-# # logging

combined

static_comb.df <- left_join(
  select(static_nocomb, cue_1, label_ci_2, log_1, Self = fitness_diff),
  select(static_comb, cue_1, label_ci_2, log_1, fTotal = fitness_diff),
  by = c("cue_1", "log_1", "label_ci_2")) %>% 
  filter(!is.na(Total) & !is.na(Self)) %>% 
  mutate(classification = ifelse(Total > Self, "Total better", "Self better"))

Error in `h()`:
! Problem with `filter()` input `..1`.
ℹ Input `..1` is `!is.na(Total) & !is.na(Self)`.
x object 'Total' not found
Backtrace:
 1. ... %>% ...
 8. base::.handleSimpleError(...)
 9. dplyr h(simpleError(msg, call))

plot

#——————————–# # invasion #——————————–#

# join invade df with label because I am lazy
invade.df2 <- invade.df %>% 
  left_join(select(ez_label, id_ci, label_ci_1 = label_ci), by = c("mut_id" = "id_ci"))

log

# get non-logged pairings
invade_nolog <- invade.df2 %>% 
  mutate(cue_1 = trimws(gsub("\\ .*", "", label_ci_1)),
         log_1 = case_when(
           str_detect(label_ci_1, "log") ~ "log",
           str_detect(label_ci_1, "log", negate = T) ~ "none")) %>% 
  filter(log_1 == "none")


invade_log <- invade.df2 %>% 
  mutate(cue_1 = trimws(gsub("\\ .*", "", label_ci_1)),
         log_1 = case_when(
           str_detect(label_ci_1, "log") ~ "log",
           str_detect(label_ci_1, "log", negate = T) ~ "none")) %>% 
  filter(log_1 == "log")

invade_log.df <- left_join(
  select(invade_nolog, cue_1, res_id, log_1, None = fitness),
  select(invade_log, cue_1, res_id, log_1, Log = fitness),
  by = c("cue_1", "res_id")) %>% 
  filter(!is.na(None) & !is.na(Log)) %>% 
  mutate(classification = ifelse(Log > None, "Logged better", "Not logged better"))

combined

invade_nocomb <- invade.df2 %>% 
  mutate(cue_1 = ifelse(
    str_detect(label_ci_1, "sum"), "I+Ig",
    trimws(gsub("+1.*|\\ log", "", label_ci_1))
  ),
  log_1 = case_when(
           str_detect(label_ci_1, "log") ~ "log",
           str_detect(label_ci_1, "log", negate = T) ~ "none"),
    comb_1 = case_when(
    str_detect(label_ci_1, "1+|sum") ~ "comb",
    str_detect(label_ci_1, "1+|sum", negate = T) ~ "none" 
  )) %>% 
  filter(comb_1 == "none")

invade_comb <- invade.df2 %>% 
  mutate(cue_1 = ifelse(
    str_detect(label_ci_1, "sum"), "I+Ig",
    trimws(gsub("+1.*|\\ log", "", label_ci_1))
  ),
  log_1 = case_when(
           str_detect(label_ci_1, "log") ~ "log",
           str_detect(label_ci_1, "log", negate = T) ~ "none"),
    comb_1 = case_when(
    str_detect(label_ci_1, "1+|sum") ~ "comb",
    str_detect(label_ci_1, "1+|sum", negate = T) ~ "none" 
  )) %>% 
  filter(comb_1 == "comb")
  
invade_comb.df <- left_join(
  select(invade_nocomb, cue_1, res_id, log_1, Self = fitness),
  select(invade_comb, cue_1, res_id, log_1, Total = fitness),
  by = c("cue_1", "log_1", "res_id")) %>% 
  filter(!is.na(Total) & !is.na(Self)) %>% 
  mutate(classification = ifelse(Total > Self, "Total better", "Self better"))

plot. nevermind just assemble in illustrator

invade_cue.pl <- ggarrange(invade_log.pl, invade_comb.pl, align = "h", ncol = 2)
ggsave(here("figures/report21/invade_cue.png"))

Saving 7 x 7 in image

#——————————–# # arrange co-infection figure #——————————–#

ggarrange(static.pl, static_cue.pl, invasion.pl, invade_cue.pl, widths = c(1.5, 1), align = "hv", labels = c("A", "", "B", ""))

ggsave(here("figures/report21/coinfection_final.png"), width = 15, height = 12)

#———————————–# # organization dual cue optimization data #————————————#

check if fitness of the dual cue is better than individual. 17 out of 25 comparisons have lower fitness lol when optimized via dual cue. Not sure what could remedy this. Smaller cue steps?

dual cue improved the fitness of the following combination. note many of these comparisons are made when time steps are different. Need to modify optimization so that time scale is the same. We should also be comparing fitness at 30 days when full resolution is achieved.

1. G_log and I_log
2. G_none, I_log
3. Ig_log and I_log
4. Ig_none, I_log
5. R_log, I_log
6. R_none, I_log
7. R_none, Ig_none

I_log is a reaccuring theme here.

dual_cue.df2 %>% 
  left_join(select(si_opt.df, id, fitness_si = fitness_20), by = c("id" = "id")) %>% 
  left_join(select(si_opt.df, id, fitness_si_b = fitness_20), by = c("id_b" = "id")) %>% 
  mutate(fitness_diff = fitness - fitness_si,
         fitness_diff_b = fitness - fitness_si_b) %>% 
  filter(fitness_diff > 0 & fitness_diff_b > 0)

plotting out some interesting heatmaps

source(here("functions/par_to_hm.R"))

# G log and I log
par_to_hm(par = c(-0.535568221, -2.5019109, 5.3812132,  1.8928449,  0.23451243),
          cue_range = seq(0,    4.77815125, length.out = 500),
          cue_range_b = seq(0,  6.77815125, length.out = 500),
          plot = T)

# Ig_log and I_log
par_to_hm(par = c(-0.253533013, -3.1078035, 7.8203134,  0.3603524,  -3.51389462),
          cue_range = seq(0,    6.77815125, length.out = 500),
          cue_range_b = seq(0,  6.77815125,length.out = 500),
          plot = T)

dual_cue.df2

# R_log, I_log
par_to_hm(par = c(5.249264033   ,8.9174937, -52.6281549 ,-21.1288712,   -13.05080084),
          cue_range = seq(6,    7,  length.out = 500),
          cue_range_b = seq(0,  6.77815125,length.out = 500),
          plot = T)

NA